Morphologies of three-dimensional shear bands in granular media 
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We present numerical results on spontaneous symmetry breaking strain localization in axisym- 
metric triaxial shear tests of granular materials. We simulated shear band formation using three- 
dimensional Distinct Element Method with spherical particles. We demonstrate that the local shear 
intensity, the angular velocity of the grains, the coordination number, and the local void ratio are 
correlated and any of them can be used to identify shear bands, however the latter two are less 
sensitive. The calculated shear band morphologies are in good agreement with those found exper- 
imentally. We show that boundary conditions play an important role. We discuss the formation 
mechanism of shear bands in the light of our observations and compare the results with experiments. 
At large strains, with enforced symmetry, we found strain hardening. 

PACS numbers: 45.70.Cc, 81.40.Jj 
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I. INTRODUCTION 

The description of the rheological properties of dry 
granular media is a key question which controls the abil- 
ity of handling (mixing, storing, transporting, etc.) of 
these particulate systems. An interesting and sometimes 
annoying feature of such materials is strain localization 
which appears almost always when a sample is subjected 
to deformation. The morphology of these narrow do- 
mains (shear bands) is far from being understood. 

Two-dimensional and boundary induced shear band 
shapes have a vast literature dating back to decades in- 
cluding both numerical and experimental studies. Three- 
dimensional studies have a major drawback in the diffi- 
culty of getting information from inside the sample. How- 
ever, in the past 20 years they gained increasing attention 
as experimental tools as Computer Tomography (CT) be- 
came available. Such experimental studies P, 0, H, 0| 
revealed complex localization patterns and shear band 
morphologies depending on the test conditions. 

In this paper we focus on triaxial tests, which in gen- 
eral are elementary tests, performed to obtain mechani- 
cal properties of soils. The most common axisymmetric 
triaxial test consists of a cylindrical specimen enclosed 
between two end platens and surrounded by a rubber 
membrane on which an external pressure is applied (see 
for example Q). The end platens are pressed against 
each other in a controlled way: Either with constant ve- 
locity (strain control) or with constant force (stress con- 
trol). The force resulting on the platens or the displace- 
ment rate of the platens is recorded as well as the volume 
change of the specimen. 

We report a numerical study of triaxial tests of co- 
hesionless granular materials based on three-dimensional 
Distinct Element Method (DEM). We show that depend- 



ing on the boundary conditions different shear band mor- 
phologies can be observed similar to experiments. We 
identify the shear bands by calculating the local shear 
intensity and show that it is correlated with the angular 
velocity of the particles and also with the local void ra- 
tio and the coordination number, which give alternative 
ways to detect shear bands and further verification of our 
results. 



II. SIMULATIONS 

We used a standard DEM with Hertz contact model 
Q and appropriate damping f(5] combined with a fric- 
tional spring-dashpot model [7, 8] . The triaxial tests were 
performed on vertical cylindrical samples (see Fig. [I} of 
diameter D = 22 mm and height H = 46 mm (i.e. the 
slenderness was H/D ss 2.1 similar to most experiments). 
Each sample consisted of 27000 spherical particles with 
the same mass density g = 7.5- 10 3 kg/m 3 . The particles 
had a Gaussian size distribution. The mean particle di- 
ameter was d = 0.9 mm. The standard deviation of the 
particle diameters was Ad = 0.025 mm. 

The normal F n and tangential F t components of the 
contact force were calculated as 



- K A3/2 _ xl/2 

= K t 5 t - 7tv t , 
N/m 3 / 2 , Kt = 10 4 



(1) 
(2) 



Ft : 

where K n = 10 6 N/m 3 / 2 , K t = 10 4 N/m, j„ 
1 N s/m 3 / 2 , and 7t = 1 N s/ m are the normal and tan- 
gential stiffness and damping coefficients, S n and S t are 
the normal and tangential displacements, and v n and v t 
are the normal and tangential relative velocities. 

The normal displacement was calculated from the rel- 
ative position, the size, and the shape of the bodies in 
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FIG. 1: (Color online) a) A granular sample was subjected to 
axial load and confining pressure, b) The rubber membrane 
surrounding the sample was simulated by overlapping spheres 
initially arranged in a triangular lattice, c) The neighboring 
spheres were interconnected with linear springs. The confin- 
ing pressure acted on the triangular facets. 



contact. The tangential displacement was calculated by 
integrating the tangential velocity in the contact plane 
during the lifetime of the contact. The Coulomb law lim- 
its the (tangential) friction force to fj,F n (where /i = 0.5 
is the used coefficient of friction). To allow for sliding 
contacts, we limited the length of the tangential displace- 
ment to fj,F n /Kt- (For a review on DEM see 0, S] and 
references therein. For more details on our implementa- 
tion see Q.) 

The translational motion of bodies is calculated with 
Verlet's leap-frog method. The rotational state is inte- 
grated in quaternion representation with Eulcr's method. 
With the above stiffness and damping coefficients, the in- 
verse of the average eigenfrequency of contacts, in both 
normal and tangential direction, is more than one order 
of magnitude larger than the used integration time step 
At = 10 -6 s. This assures that the noise level induced 
by numerical errors and grain elasticity is kept low. 

The initial configuration was generated by randomly 
placing the grains in a tall solid cylinder having height 
h = 3H and width D. The maximum allowed initial 
grain overlap was 1%. The upper platen and the parti- 
cles were given downward velocities v oc z/h depending 
on their vertical position z measured from the fixed bot- 
tom platen. These conditions lead to an almost simulta- 
neous first contact of the bodies. The system was stabi- 
lized with a force applied on the upper platen. This was 
switched on when the inner pressure of the sample could 
compensate it. In preparation (and later in the tests) 
we used zero gravity. The deposition method described 
above is known to produce a homogeneous system. Dur- 
ing the preparation phase, the friction of the particles 
was switched off to allow for generation of dense sam- 
ples. (For a review on sphere packings see [T(J.) 

The solid cylinder used in preparation was replaced 
by an elastic membrane in the tests. The elastic mem- 
brane was modeled with overlapping spheres having equal 
diameter d m — 1 mm and equal mass density g m — 



100 kg/m 3 , and initially forming a triangular lattice on 
the external surface of the cylinder. The "membrane 
nodes" could not rotate and were interconnected with 
linear springs having an elongation equal to the relative 
distance of the nodes (initially 0.5 mm). The stiffness 
of the springs k s = 0.5 N/m was chosen such that the 
particles could not escape by passing through the mem- 
brane. Additionally a homogeneous confining pressure 
<r c = 500 N/m 2 was applied on the membrane. This was 
simulated by calculating the forces acting on the trian- 
gular facets formed by connected membrane nodes. 

A similar model was used by Tsunekawa and Iwashita 
PH who applied the confining pressure directly on the 
external particles in a very similar way. However, their 
approach requires the computationally expensive identifi- 
cation of external particles and a Delaunay triangulation. 
Sakaguchi and Muhlhaus [l2j used a similar membrane 
model to ours but without an explicit confining pressure, 
relying only on the stiffness of the springs. 

The bottom platen was fixed during both preparation 
and test phases. The upper platen could not tilt in 
preparation phase, but in certain tests it could freely 
tilt along any horizontal axis with rotational inertia 
/ = 10 -7 kg m 2 . During the tests, the samples were 
compressed by moving the upper platen in vertical direc- 
tion downwards with a constant velocity (strain control) . 
Starting from the same initial condition, four different 
runs - denoted by (A), (B), (C), and (D) - were exe- 
cuted. Two different compression velocities were used: 
A base value u\ = 10 mm/ s (in tests (A) and (B)) and 
a two times larger value u-i = 1u\ (in tests (C) and (D)). 
Tilting of the upper platen was enabled in tests (A) and 
(C) and disabled in tests (B) and (D). 



III. RESULTS 

A. Local shear intensity 

We define the local shear intensity S by generalizing its 
two-dimensional definition given by Daudon et al. [1 31 ] - 
First, the regular triangulation [3] of the particle system 
is calculated [15j . The displacements of the particles (rel- 
ative to a previous state) are known from the DEM sim- 
ulation. We extend the displacement field to the whole 
volume of the sample with a linear interpolation over the 
tctrahedra of the regular triangulation. For each particle 
we identify the incident triangulation cells (tetrahedra) 
which define a discrete particle neighborhood. The par- 
ticles at the sample's boundary, having infinite incident 
cells, are skipped (i.e. no local shear intensity is defined 
for them). The discrete neighborhood of a particle is a 
surrounding polyhedron with first-neighbor particles at 
the corners. 

We define the deformation gradient tensor with the 
partial (space) derivatives diUj of the displacement vector 
u. In a neighborhood fl of volume V, the components of 
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the mean deformation gradient tensor are calculated as 



1 



u v) = y I diUjdV. 



(3) 



Using the Gauss-Ostrogradski theorem the volume inte- 
gral can be transformed into a closed surface integral over 
the boundary dfi of f2, leading to 



( u ij) = 77 / n i u j dS, 
V Jan 



(4) 



where n is the exterior normal along the boundary. We 
calculate the local deformation gradient tensor applying 
the above formula to discrete neighborhoods and using 
the linear interpolation of the particle displacements. In 
this case the integral can be reduced to a summation over 
triangular facets. 

The symmetric part of the local deformation gradient 
tensor is a macroscopic strain tensor derived from parti- 
cle displacements. Using the eigenvalues Ek of this macro- 
scopic strain tensor, we define the local shear intensity as 



S 



max 

k 



(5) 



We note that we disregard the elastic deformation and 
rotation of the grains, since we are interested in the iden- 
tification of the shear bands, which is strongly linked to 
geometric effects. However, for constitutive models, a 
more complete treatment of the strain would be needed 



B. Shear band morphologies 

Taking cross sections of the sheared samples and col- 
oring the grains according to the local shear intensity 
S, we could identify shear bands (see Fig. [5]) and com- 
pared them with experiments. In experiments the CT 
scans show the volume fraction difference between the 
bulk and the shear band. In the next section we justify 
the comparison of the volume fraction and local shear 
intensity. 

Our simulations are run at zero gravity and low confin- 
ing pressure in very similar conditions to the experiments 
of Batiste et al. Q. The shear band patterns found in 
their experiments and our simulations are also very close 
to those found in experiments under normal gravity and 
high confining pressure by Desrues et al. [2|, who also 
studied the case of a tilting upper platen. We compared 
our results to both experiments. 

We found that the absence of enforced axisymmetry 
leads to spontaneous symmetry breaking. When tilting 
of the upper platen is enabled internal instabilities can 
develop into a localized deformation along a failure plane 
(see panel (c) of Fig. [5]). Non-tilting platens act as a 
stabilizing factor leading to an axisymmetric hourglass 
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FIG. 2: (Color online) Cross sections (a, c, e, f, g) of sample 
(C) and cross sections (b, d, h, i, j) of sample (D) shown 
at 10% axial strain. Panels (k, 1, m) present CT scans 0] 
(F2075). The vertical cross sections (a-d) were taken at the 
middle of the sample. The horizontal cross sections were taken 
at different heights: close to the top (e, h, k), at the middle 
(g, j, m), and in between (f, i, 1). The red color encodes the 
local shear intensity on panels (c-j) and the local void ratio 
on panels (k-m). Panels (a, b) show the velocity field. 



shaped shear band with two conical surfaces and com- 
plex localization patterns around them (see panel (d) of 
Fig. [2]). This is in full accordance with the experimen- 
tal results of Desrues et al. 0]. For the non-tilting case 
Tsunekawa and Iwashita [ll[ found in DEM simulations 
similar localization patterns, however they have not in- 
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vestigated the tilting case. 

In the non-tilting case, the horizontal cross sections (h, 
i, j) shown on Fig. [5] can be compared with the exper- 
imental results of Batiste et al. j3j]. They reported the 
same type of shear band morphologies for these boundary 
conditions. The panels (k, 1, m) of Fig. [2]show CT scans 
from their triaxial shear tests executed in micro-gravity 
aboard a NASA Space Shuttle. Our simulations used 
similar setup and similar confining pressure. Good agree- 
ment of shear band shapes (including their non-trivial 
structure) can be recognized in spite of the rather lim- 
ited number of grains in our simulations. (Note that the 
details reproduce better in the color version of the figure.) 

In the case of tilting upper platen the shear bands are 
not totally plane - as can be seen on horizontal cross sec- 
tions taken close to the platens - but follow the curvature 
of the boundary. The same was found experimentally by 
Desrues et al. Q ■ In the vertical cross sections shown on 
the panels (c) and (d) of Fig. [2 the found shear bands 
are in good agreement with changes in the velocity field 
shown on panels (a) and (b). This justifies the shear band 
identification method based on the local shear intensity. 
We have also investigated alternative methods. 



C. Alternative methods of shear band identification 

It is widely known that dense granular materials dilate 
during shear. In some experiments (e.g. experiments 
based on CT 0, H[) the local void ratio is used to iden- 
tify the shear bands. To confirm the presented shear 
band identification method and the found shear band 
morphologies, we have investigated the correlation be- 
tween the local void ratio v and the local shear intensity 
S. The void ratio was measured using the regular trian- 
gulation [3 EH of the spherical particles. The volume of 
the regular Voronoi cells V c and the volume of the grains 
V g define the local void ratio v = (V c — V g )/V g . 

In numerical simulations (for spheres of nearly equal 
size), a good alternative to the local void ratio v is the 
coordination number Z (defined by the number of con- 
tacts), which decreases as v increases. Its main advantage 
is that it can be defined exactly and calculated fast, how- 
ever, if the size distribution is wide a non-trivial particle 
size scaling has to be taken into account. 

The existence of particle rotations in shear bands is 
known to experimentalists for a long time (see e.g. [I?) )• 
It was also evidenced in simulations by Herrmann et al. 
[III]. In our simulations, we have also measured for each 
grain the absolute value of the angular velocity R and 
tested its correlation with the local shear intensity S. 

All the quantities mentioned above (the local void ra- 
tio v, the coordination number Z, the angular velocity 
R, and the local shear intensity S) are defined for each 
particle. We checked their correlation with a histogram 
technique using the values calculated for different par- 
ticles as different statistical samples. The v, Z, and R 
values were averaged for each (logarithmic) histogram bin 
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FIG. 3: Correlation of the local shear intensity S with the 
local void ratio [y, o), the coordination number (Z, x), and 
the angular velocity of the grains (R, +). X is one of v, Z, and 
R. F denotes a linear transformation different for each data 
set. All quantities are scaled by average values (X av , S a v)- 
The data is collected from four samples at 10% axial strain. 
See text for more details. 

of S. We also calculated the total averages v av , Z av , R av , 
and S av . On the quantities £ = ln(X / X av ) (where X is 
one of v, Z, and R) we applied different linear transfor- 
mations F(£) = a(£ — £o) (shift and scaling) to achieve 
data collapse of F(\n(X/X av )) as function of \n(S/S av ) 
(see Fig. Ej. 

The scaling term a of F shows the sensitivity of the R, 
Z, and v quantities with respect to the local shear inten- 
sity S. We found a = 1 for the angular velocity, a = —9 
for the coordination number, and a = 27 for the local 
void ratio. (Note that Z decreases as S increases!) The 
fluctuations were proportional to \f\u\. Regarding shear 
band identification, this means that the angular veloc- 
ity is essentially equivalent with the local shear intensity. 
However, the coordination number and the local void ra- 
tio are less sensitive, and they exhibit large fluctuations 
due to random packing and random rearrangements. For 
this reason they need more spatial and/or temporal av- 
eraging to achieve the same accuracy. 



D. Stress-strain relation 

In order to compare to most common experimental re- 
sults, we measured the stress a on the upper platen, and 
calculated the stress ratio a/ao, where oo denotes the 
initial stress. As the axial strain increases, the response 
of the granular sample (the stress ratio) increases until 
it reaches a peak value, then it decreases (see Fig. [4j. 
According to Fig. [H up to 15% axial strain there is no 
significant difference in stress-strain relation measured in 
different simulation runs, indifferent to strain rate and 
tilting of the upper platen. (We have not tested depen- 
dence on material parameters and confining pressure.) 

The presented strain softening effect is a basic obser- 
vation of triaxial shear tests of dense granular specimens 
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FIG. 4: (Color online) Stress-strain relation. The stress ra- 
tio a/ (To (where <jq denotes the initial stress) measured on 
the upper platen is shown as function of the axial strain, for 
different simulation runs. (See inset for low strains.) For the 
lower two curves (a, c) tilting of the upper platen was enabled, 
and for the upper two (b, d) it was disabled. For (c, d) the 
samples were compressed two times faster than for (a, b). 




FIG. 5: Experimental stress-strain relation from a triaxial 
compression test exhibiting strain softening and development 
of shear plane. (Courtesy of Poul V. Lade, Reprinted from 
with permission from Elsevier.) 

(see for example [ll| and Fig. [5]). Any local deforma- 
tion due to shear is followed by a dilatation resulting in 
a decrease of the force bearing capacity of the material 
which further intensifies the deformation leading to fail- 
ure. In our simulations we observed pronounced shear 
bands around 10% axial strain, which is in good agree- 
ment with the experimental results. 

After about 15% axial strain the different bound- 
ary conditions and shear band shapes result in different 
stress-strain curves. In the non-tilting case, the geometry 
and the hourglass-shaped shear band forces the particles 
to enter and leave the failure zone. In this case, a stable 



slipping mode cannot be formed. As the test sample is 
further compressed, it opposes more and more firmly to 
compression (see curves (b) and (d) on Fig. 0]), result- 
ing in increasing stress ratio (i.e. strain hardening). In 
the tilting case, the upper part of the sample moves as 
a single block. The formed planar shear band allows for 
a stable slipping mode with nearly constant stress until 
boundary effects come into play (see curves (a) and (c) 
on Fig. 111). 



IV. CONCLUSIONS 

Triaxial shear test simulations based on DEM were exe- 
cuted and different shear band morphologies known from 
experiments were reproduced [2lj . We have shown that 
in triaxial shear tests symmetry breaking strain localiza- 
tion can develop spontaneously if the axial symmetry is 
not enforced by non-tilting platens. To our knowledge 
it is the first time that such symmetry breaking strain 
localization was reproduced in DEM simulations. 

We generalized the shear intensity definition of Daudon 
et al. [HI to three-dimensions and used it to identify 
shear bands. To be able to compare to experiments, we 
have also tested alternative methods of shear band iden- 
tification. We found strong correlation of the local shear 
intensity with the angular velocity of the grains, the co- 
ordination number, and the local void ratio. This result 
justifies our method and proves once more the known ex- 
perimental and numerical findings that shear bands are 
characterized by dilation and rotation of the grains. Re- 
garding shear band identification, the coordination num- 
ber and the local void ratio are found to be less sensitive 
than the local shear intensity and the angular velocity of 
the grains. 

We have also measured the stress-strain relation of the 
compressed samples. Strain softening was identified in 
good agreement with experimental results. We have also 
found a strain hardening effect at large strains in the non- 
tilting case and explained it in terms of geometry and 
shear band morphology. However, this might be only 
valid for the tested material parameters and confining 
pressure. We have no knowledge of experiments focusing 
on this particular question. In general, the agreement of 
our results with the experimental results is very good, 
even if the system size (number of particles) in our sim- 
ulations is much smaller than in experiments. 
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